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ABSTRACT 

Although much research has been devoted to the problem of 
restoring Poissonian images, namely in the fields of medical 
and astronomical imaging, applying the state of the art regu- 
larizes (such as those based on wavelets or total variation) to 
this class of images is still an open research front. This pa- 
per proposes a new image deconvolution approach for images 
with Poisson statistical models, with the following building 
blocks: (a) a standard regularization/MAP criterion, combin- 
ing the Poisson log-likelihood with a regularizer (log-prior) 
is adopted; (b) the resulting optimization problem (which is 
difficult, since it involves a non-quadratic and non-separable 
term plus a non-smooth term) is transformed into an equiva- 
lent constrained problem, via a variable splitting procedure; 
(c) this constrained problem is addressed using an augmented 
Lagrangian framework. The effectiveness of the resulting al- 
gorithm is illustrated in comparison with current state-of-the- 
art methods. 

1. INTRODUCTION 
1.1. Poissonian Images 

Image restoration is one of the earliest and most classical in- 
verse problems in imaging, dating back to the 1960's. Much 
of the work in this field has been devoted to developing reg- 
ularizers (priors or image models, in a Bayesian perspective) 
to deal with the ill-conditioning or ill-posedness of the obser- 
vation operator, and to devising efficient algorithms to solve 
the resulting optimization problems. 

A large fraction of the work on image restoration assumes 
that the observation operator is linear (often the convolution 
with some blur point spread function) and the presence of 
additive Gaussian noise. For this scenario, recent work has 
lead to a set of state-of-the-art restoration methods, which in- 
volve non-smooth convex regularizes {e.g., total-variation, £i 
norm of frame coefficients) and efficient special-purpose al- 
gorithms (see J2], 0, lfl2l . Q3Q , and references therein). 

The algorithms developed for the linear/Gaussian obser- 
vation model cannot be directly applied to other statistical 
(e.g., Poisson or Gamma) observation models. The Poisson 



case is well studied and highly relevant in fields such as as- 
tronomical [18], biomedical ||8), JTT |, and photographic imag- 
ing iTOl . A very recent overview of deconvolution methods 
for Poissonian images can be found in [9 |, where a state-of- 
the-art algorithm is also introduced. 

Although our approach can be applied to other regular- 
izes, we focus here on total-variation (TV), well-known for 
its discontinuity preserving ability [3], [16|. The combina- 
tion of TV regularization with the log-likelihood resulting 
from the Poissonian observations of a convolved image, leads 
to an objective function with a non-quadratic non-separable 
term (the log-likelihood) plus a non-smooth term (TV). This 
objective function poses the following difficulties to the cur- 
rent state-of-the-art algorithms: (a) the Poisson log-likelihood 
term doesn't have a Lipschitz-continuous gradient, which is a 
necessary condition for the applicability of algorithms of the 
forward-backward splitting (FBS) class Q, (9); (b) the pres- 
ence of a convolution in the observation model precludes the 
direct application of the Douglas-Rachford splitting methods 
described in [6|. Moreover, if an FBS algorithm is applied 
(ignoring that the convergence conditions are not met), it is 
known to be slow, specially when the observation operator is 
ill-conditioned, a fact which has stimulated recent research 
aimed at obtaining faster methods JT], 0, I2TI . 

In this paper, we propose a new approach to tackle the 
optimization problem referred to in the previous paragraph. 
Firstly, the original optimization problem is transformed into 
an equivalent constrained one, via a variable splitting pro- 
cedure. Secondly, this constrained problem is addressed us- 
ing an algorithm developed within the augmented Lagrangian 
framework, for which convergence is guaranteed. The effec- 
tiveness of the resulting algorithm is illustrated in comparison 
with current state-of-the-art alternatives (9), iTOl . 0. 



2. AUGMENTED LAGRANGIAN 

In this section, we briefly review the augmented Lagrangian 
framework, a key building block of our approach. Consider a 



convex optimization problem with linear equality constraints 
min E(v) 

vef (1) 
s.t. Av = b, 

where b G R p and A G R pxd . The so-called augmented 
Lagrangian function for this problem is defined as 

C A (v, V, M) = ^(v) + »7 T (Av - b) + | || Av - (2) 

where r\ G W is a vector of Lagrange multipliers and /i > 
is called the AL penalty parameter |15|. The AL algorithm 
iterates between minimizing £a(v, with respect to v, 
while keeping r\ fixed, and updating rj. 

Algorithm AL 

1. Set k = 0, choose /i > 0, vo, and T] . 

2. repeat 

3. v fc+ i G argmin v £ j4 (v, rj fc ,/i) 
4 - »7fc+i <- *7fc + M(Av fc+ i - b) 

5. + 1 

6. until stopping criterion is satisfied. 

It is possible (in some cases recommended) to update the 
value of /i at each iteration 1151 . Notice, however, that it is 
not necessary to take /i to infinity to guarantee convergence 
to the solution of the constrained problem dTJ. In this paper, 
we will consider only the case of fixed /i. 

After a straightforward manipulation, the terms added to 
E(v) in C^(v,T] k , fi) (see (O) can be written as a single 
quadratic term, leading to the following alternative form for 
the AL algorithm: 

Algorithm AL (version 2) 

1. Set k — 0, choose /i > 0, vo, and do. 

2. repeat 

3. v fc+ i e argmin v £;(v) + £||Av - d fc ||| 

4. d fc+ i <- d fc - (Av fe+1 - b) 

5. k^k + 1 

6. until stopping criterion is satisfied. 

This form of the AL algorithm makes clear its equivalence 
with the recently introduced Bregman iterative method l22l . 

3. PROBLEM FORMULATION 

Let y = (j/i, y n ) G Nq denote an n-elements observed 
image or signal of counts, assumed to be a sample of a random 
image Y = (Yi, ...,Y n ) G Nq composed of n independent 
Poisson variables 

P[Y = y|A]=n^-r-, C3) 

i=i yt ' 



where A = (Ai, A n ) G IR™ is the underlying mean signal, 
assumed to be a blurred version of an unknown x, i.e., 

A = K x, (4) 

where K is the matrix representation of the blur operator, 
which is herein assumed to be a convolution. When dealing 
with images, we adopt the usual vector notation obtained by 
stacking the pixels into an n-vector using, e.g., lexicographic 
order. Combining (fj) and we can write 

n 

logP[Y - y|x] = Y^V* M(Kx)0 - (Kx) f - log^!) 

i=i 

where (Kx)i denotes the i-th component of K x |8|, |18|. 

Under the regularization or the Bayesian maximum a pos- 
teriori (MAP) criterion, the original image x is inferred by 
solving a minimization problem with the form 

min L(x) (5) 
s.t. x > 0. (6) 

The objective L(x) is the penalized negative log-likelihood, 

L(x) = -lo g P[Y = y|x]+T0(x), (7) 
= ^(Kx),- !/l lo 6 ((Kx) 1 )+r ( t(x), (8) 

i=l 

where <f> '■ R n —* R is the penalty/regularizer (negative of 
the log-prior, from the Bayesian perspective), and r G M+ is 
the regularization parameter. Notice that the non-negativity 
constraint on x guarantees that A = K x is also non-negative, 
if all the entries in K are non-negative (as is the case in most 
convolution kernels modeling a variety of blur mechanisms). 
In this work, we adopt the TV regularizer 0, fl6l . i.e., 

" / 

0(x) = TV(x) - V ( A ^ x ) 2 + ( A ^ x ) 2 ' ^ 

s=l 

where ( x and A^x) denote the horizontal and vertical first 
order differences at pixel s G {1, . . . , n}, respectively. 

Each term (K x)^ — log ((K x)^) of ([S]), corresponding 
to the negative log-likelihood, is convex, thus so is their sum. 
If the space of constant images {x = a(l, 1, 1), a G M}, 
for which TV is zero, does not belong to the null space of K, 
and the counts y n ) are all non-zero, then the objective 

function L is coercive and strictly convex thus possessing a 
unique minimizer Q . 

4. PROPOSED APPROACH 
4.1. Variable Splitting 

The core of our approach consists in rewriting the optimiza- 
tion problem defined by ©-(O as the following equivalent 



constrained problem: 



mm 

x,z,u 



S. t. 



^2(zi -yi logs,-) +T<f>(\l) (10) 

i=l 

Kx = z (11) 

x = u. (12) 



Notice that we have dropped the non-negativity constraint ©; 
this constraint could be applied to either x, z, or u (as long 
as all elements of K are non-negative). However, as shown 
below, if applied to z, this constraint will be automatically 
satisfied during the execution of the algorithm, thus can be 
dropped. Notice that this problem can be written compactly 
in the form (Q]), using the translation table 



x 
z 
u 



b = 0, A = 



K I 
I I 



(13) 



and with 

n 

E(v) = E(x,z,u) = J2( z i-Vi logz,) +t0(u). (14) 



4.2. Applying the AL Algorithm 

The application of Step 3 of the AL (version 2) algorithm 
to the problem just described requires the solution of a joint 
minimization with respect to x, z, and u, which is still a 
non-trivial problem. Observing that each partial minimiza- 
tion (e.g., with respect to x, while keeping z and u fixed) is 
computationally treatable suggests that this joint minimiza- 
tion can be addressed using a non-linear block Gauss-Seidel 
(NLBGS) iterative scheme. Of course, this raises the question 
of wether such a scheme converges, and of how much compu- 
tational effort {i.e., iterations) should be spent in solving this 
minimization in each step of the AL algorithm. Experimental 
evidence (e.g. [14|) suggests that good results are obtained 
by running just one NLBGS step in each step of the AL algo- 
rithm. In fact, it has been shown that the AL algorithm with 
a single NLBGS step per iteration does converge iflOll . ifTTl . 
Remarkably, the only condition required is that the objective 
function be proper and convex. 

Finally, applying AL (version 2), with a single NLBGS 
step per iteration, to the constrained problem presented in the 
previous subsection leads to our proposed algorithm, termed 
PIDAL (Poisson image deconvolution by AL). The algorithm 
is presented in Fig.Q] 

The minimization with respect to z (line 5) is given by 



_ 1 = (K T K + l)" 1 (K T x' + x"). 



(15) 



We are assuming that K models a convolution, thus it is a 
block Toeplitz or block circulant matrix and ( TT3T > can be im- 
plemented in 0(n log n) operations, using the FFT algorithm. 



Algorithm Poisson Image Deconvolution by AL (PIDAL) 

1. Choose xo, zo, uo, d^, dp 2 ', /x, and r. Set k := 0. 

2. repeat 



3. 
4. 
5. 

6. 
7. 

8. 
9. 

10. 

11. 
12. 



Zfc 



x" = u fc + df 

x fc+ i := argmin ||Kx 

X 

z' = Kx fe+1 - d< 1} 



x'Hl + ||x - 



x"||l 



z fc+ i := arg mm 



Vi logZ t + — || Z - Z || 2 



U = X fc+ i - vi fc 

■ 1,1 ',,2 

u fc+ i := arg mm - u - u 

x Z 

dj& :=dl l3 -(Kx fe+1 -z fc+1 ) 



|| 2 + (r/ M )^(u). 



,(2) 



- fc+1 :=d^ 2) - (xfc+i -u fc+ i) 
k := k + 1 



13. until some stopping criterion is satisfied. 



Fig. 1. The PIDAL algorithm. 



Step 7 is separable and has closed form: for each Zi, it 
amounts to computing the non-negative root of the second 
order polynomial fizf + (1 — /i z'^Zi — y,, given by 



1 + ((M- 



1) 



4 n y i 



,1/2 



/(2M). 
(16) 

Notice that this is always a non-negative value, thus justifying 
the statement made above that the constraint z > is auto- 
matically satisfied by the algorithm. 

The minimization with respect to u (line 9) is, by defini- 
tion, the Moreau proximity mapping : K™ — » K™ of the 
regularizer T<j> Q- In this paper, the adopted regularizer is the 
TV norm (O, thus u^+i is obtained by applying TV-based de- 
noising to u'. To implement this denoising operation, we use 
Chambolle's well-known algorithm [3], although other fast 
methods are also available [20|. 

Notice how the variable splitting, followed by the aug- 
mented Lagragian approach, converted a difficult problem (|5]l- 
(O, involving a non-quadratic and non-separable term plus a 
(non-smooth) TV regularizer, into a sequence of three sim- 
pler problems: (a) quadratic problem with a linear solution 
(line 5); (b) a separable problem with closed-form solution 
(line 7); (c) a TV-based denoising problem (line 9), for which 
efficient algorithms exist. 

5. EXPERIMENTS 

We now report experiments where PIDAL is compared with 
two state-of-the-art methods f9), |[T3l . All the experiments 
use synthetic data produced according to ©-(lU, where x 
is the Cameraman image and K represents a uniform blur. 
In Experiment 1 (following |[T3l ). the blur is 9 x 9, and the 



Table 1. Mean absolute errors obtained by PIDAL and the 
algorithm from [9] (average over 10 runs). 



max intensity 


5 


30 


100 


255 


PIDAL 


0.37 


1.34 


3.99 


8.65 


Algorithm from [9 1 


0.44 


1.44 


4.69 


10.40 



original image is scaled to a maximum value of 17600; this 
is a high SNR situation. In Experiment 2 (following (9)) 
the blur is 7 x 7, and the maximum value of x belongs to 
{5, 30, 100, 255}; this represents low SNR situations. 

Parameter fi of the PIDAL algorithm affects its conver- 
gence speed, but its adaptive choice is a topic beyond the 
scope of this paper. In all the experiments, we use [i = r/50, 
found to be a good rule of thumb. PIDAL is initialized with 
x = y, z = K x , u = x , = 0, and d{, 2) = 0. 

In Experiment 1, the regularization parameter r was set to 
6 x 10~ 4 ; since our goal is to propose a new algorithm, not a 
new deconvolution criterion, we didn't spend time fine tuning 
r or using methods to adaptively estimate it from the data. 
Since the method in [ 1 3 1 includes a set of adjustable param- 
eters which need to be hand tuned, the comparison remains 
fair. The improvement in SNR (ISNR) obtained by PIDAL 
was 6.96dB (average over 10 runs), better than the 6.61dB 
reported in [13|. This result is more remarkable if we no- 
tice that the TV regularizer is considerably simpler than the 
locally adaptive approximation techniques used in ||T3l . 

For Experiment 2, we downloaded the code available at 
www . greyc . ensicaen . f r/~f dupe/ . Although the regu- 
larizer is not the same, we used the same values of r found 
in that code; if anything, this constitutes a disadvantage for 
PIDAL. Following [9|, the accuracy of an image estimate x 
is assessed by the mean absolute error MAE = ||x — xlli/n. 
Table Q] shows the MAE values achieved by PIDAL and the 
algorithm of |9 1, for the several values of the maximum origi- 
nal image intensity, showing that PIDAL always yields lower 
MAE. In our experiments, each run of the algorithm from [9 1 
takes roughly 10 times longer than PIDAL. 

6. CONCLUDING REMARKS 

We have proposed an approach to TV deconvolution of Pois- 
sonian images, by exploiting a variable splitting procedure 
and augmented Lagrangian optimization. In the experiments 
reported in the paper, the proposed algorithm exhibited state- 
of-the-art performance. We are currently working on extend- 
ing our methods to other regularizers, such as those based on 
frame-based sparse representations. 
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